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Abstract 

The self-diffusion coefficient of a granular gas in the homogeneous cooling state is analyzed 
near the shearing instability. Using mode-coupling theory, it is shown that the coefficient diverges 
logarithmically as the instability is approached, due to the coupling of the diffusion process with 
the shear modes. The divergent behavior, which is peculiar of granular gases and disappears in 
the elastic limit, does not depend on any other transport coefficient. The theoretical prediction is 
confirmed by molecular dynamics simulation results for two-dimensional systems. 
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Granular gases are fluidized systems composed by particles colliding inelastically. In spite 
of some apparent similarities, they behave very differently from molecular fluids, and exhibit 


many interesting and peculiar phenomena [1]. This includes the spontaneous development of 
strong density and temperature inhomogeneities, spontaneous symmetry breaking, pattern 
formation, and segregation in systems composed of different types of particles, to mention 
a few examples Here, another interesting feature of granular gases is addressed, namely 
the divergence of the Navier-Stokes self-diffusion coefficient when the shearing instability is 
approached, indicating the existence of anomalous diffusion. This effect is quite different 
from the divergence of the self-diffusion coefficient in an infinite molecular two-dimensional 
system, as a consequence of the algebraical long time tails of the velocity autocorrelation 
function (VACF) 3, 4]. The divergence being discussed here occurs in rather small systems 
and it is due to the singular behaviour associated with the presence of a hydrodynamic 
instability. 

In this paper, a granular gas is modeled as a system of N identical smooth inelastic 
hard spheres or disks of mass m and diameter a, enclosed in a volume V, and with the 
collisions characterized by a constant coefficient of normal restitution a. As a consequence 
of the energy dissipation in collisions, granular gases are always out of equilibrium. For 
isolated systems, there is a time-dependent homogeneous cooling state (HCS), in which the 


granular temperature T(t) decreases in time following Haff’s law Q, d t T(t) = -C (t)T(t), 
C(t) oc T(t) 1 / 2 being the cooling rate. This is the reference state used to derive hydrodynamic 
equations for granular gases [6H8] . The predictions from these macroscopic equations are in 
good agreement with simulation and experimental results for simple situations and dilute 


or moderately dense systems { 2 }]. One important result from hydrodynamics, confirmed by 
numerical simulations, is that the HCS is unstable with regards to spatial perturbations with 
wavelength larger than some critical value, that depends on the inelasticity of the system 


JQ. 


The origin of the instability is in the fluctuations of the shear mode that lead, through 


nonlinear hydrodynamic contributions, to the development of density inhomogeneities 
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Self-diffusion is the prototype transport process, and the associated diffusion equation is 
the prototype hydrodynamic equation for a macroscopic description of the process. It has 


been extensively investigated 
granular gases in the HCS 


Dot 
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l in molecular and granular fluids and, in particular, in 


17]. In these studies, the hydrodynamic diffusion equation 
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has been derived with a prediction for the self-diffusion coefficient. It has been shown 15] 
that the coefficient, if it exists, is given by a Green-Kubo like expression. The analysis is 


simplified by a change from the original time scale t to a new one s defined by 

t 


18 


19] 


u 0 s = In 


to ’ 


( 1 ) 


where to and ujo are two arbitrary constants. In this time scale, the original time-dependent 
HCS is exactly mapped on a steady state, whose granular temperature is 
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( 2 ) 


Here, ( is the time-independent cooling rate, 


c = c(i) 


(3) 


2 v 0 (t) ’ 

with i>o (i) = (2T/m) 1//2 . The diffusion coefficient D in the steady representation of the HCS 
is identified from the diffusion equation for the number density n of tagged particles, 

dm 


ds 


= —D\7n. 


At a formal level, D is given by the Green-Kubo expression [1 


D = 


dsC(s ), 


(4) 


(5) 


that involves the VACF in the modified dynamics C(s) for a tagged particle, defined as 

C(s) = • ™(0)), (6) 

where d (2 or 3) is the dimension of the system, w(s) is the velocity of the particle in the 
modified dynamics at time s, and the angular brackets denote average with the steady HCS 
distribution. Spatial perturbations of the transversal velocity with wavevector smaller than 

1/2 

(7) 


k i = 


VstC 


V 


lead to the shearing instability of the HCS mentioned above. In Eq. (!7j) . v st = ( 2 T st /m ) 
and rj is related to the shear viscosity r/ through 

1/2 


y/2 


T. 


V = 


St 


V 


mn \ T(t) 


( 8 ) 
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It follows that systems whose linear size L is larger than L c = 2ir/k± will exhibit the shearing 
instability. For this reason, the previous studies of self-diffusion in the HCS we are aware 
of 
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restrict themselves to situations in which L is much smaller that L r . On the 


other hand, here the focus is on the peculiarities of the process when the shearing instability 
is approached. The VACF can be split into a part associated to the fast relaxation of the 
kinetic modes, followed by the contribution of the slow relaxation of the hydrodynamic type, 


C(s) — Ckin(s) + C , hyd('S). 


(9) 


In Ref. 


20], an expression for the hydrodynamic component was derived by using mode¬ 


coupling theory. The result reads 

(d - 1 )T st 
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hyd 
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mnV d 
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The summation extends over values of k compatible with the employed periodic boundary 
conditions, and in the interval k m < k < Ism , where k m = 2i t/L and Cv is of the order 
of 2n times the inverse of the mean free path. The presence of the factor k 2 /(k 2 — k\) 
is peculiar of granular gases, becoming unity in the elastic limit, in which k± = 0. It is 
associated with the existence of the shearing instability. The above expression implies that 
there is a time window over which the VACF exhibits a power-law decay on the s scale 


20]. Besides, 


( 11 ) 


having some similarities with the long time tails occurring in molecular fluids 
for s So = L 2 /4:7r 2 (rj + D ), the decay of the VACF has the exponential form 23] 

^ . . 2(d-l)f st e- s N k ™- k V +k ™ S \ 

Chl ' d(s) ~ —^nV 1 - [L/L c ) 2 ’ 

showing that the amplitude of the decay diverges as the critical length L c is approached. In 
the following, the implications of Eq. (HU|) on the behaviour of the self-diffusion coefficient 
near the instability will be discussed. Let us insist on the fact that the analysis here will 
be for finite size L < L c , contrary to the analysis leading to the usual long time tails in 
molecular gases, requiring an infinite system ji ]. Long-time tails in freely cooling granular 
gases under conditions for which the HCS is unstable have been investigated in ref. 24], 


while in ref. 25| the case of a driven granular fluid is considered. The decomposition in Eq. 
(|9]) allows us to write 

D = D kin + -Dhyd, (12) 
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where Z\ in is the contribution from the kinetic part of C(s) and 


-Dhyd — 
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mnVd 


E 


(0 


k 2 


(. k 2 -k\) ?j(k 2 -k 2 ± ) + k 2 D 
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Upon deriving this expression, it has been taken into account that all the terms on the right 
hand side of Eq. (fTUD decay exponentially in time if the system is in the parameter region in 
which the HCS is stable. Consider now the proximity of the shearing instability, i.e. that L 
is close to (and below) L c or, equivalently, that k m is close to (and above) k±. As long as the 
viscosity does not diverge as approaching the instability, it is rj(k 2 — k\) <C k 2 D. Moreover, 
it will be shown below (see Eq. (TTB1) ) that D hy d diverges as the instability is approached, 
while ZA in is expected to remain hnite. In this way, Eqs. dT^]) and (1THT) in the vicinity of the 
instability reduce to 


D ~ D hyd - 


(d - 1 )T st 


v(') 


m7iVdD hyd - ™ _ ™_l 


(14) 


h 2 — h 2 

IV IV | 

The summation over k can be carried out by using the continuous limit. In the following, a 
two-dimensional system (d = 2) will be considered. Then, one easily gets 

1/2 


D ~ .Dhyd ^ 


T st 
8Trmn / 
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(15) 


indicating a logarithmic divergent behavior of the self-diffusion coefficient as the shearing 
instability of the HCS is approached. Note that this behaviour differs from the one obtained 
by considering only the long time exponential contribution given by Eq. m ■ 

Molecular dynamics (MD) simulations of a system of inelastic hard disks in a square box 
with periodic boundary conditions have been performed. To compare the simulation results 
with the above theoretical predictions, a quite precise value of the critical size L c is needed 
for each set of parameters. In principle, L c is determined by Eq. ©, and the expressions 
for the shear viscosity and the cooling rate as obtained by using the Enskog theory can 
be used ( 21 1. Nevertheless, those expressions correspond to bare quantities, i.e. neglecting 
hydrodynamic contributions, that can be quite relevant in the vicinity of the instability. 
Then, the critical sizes of the systems to be used in the analysis of the instability region 
have been identified from the simulations. A convenient method for this is based on the 
increase of the average kinetic energy of the system near the instability. Using fluctuating 
hydrodynamics, it has been shown that the average of the kinetic energy per particle (e) in 
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FIG. 1: (Color online) Critical size L c for the shearing instability as a function of the number 
density for two values of the coefficient of normal restitution a. The symbols are simulation results 
using Eq. (1161) . while the lines are the theoretical predictions given by Eq. ([7]). using the Enskog 
values for the cooling rate and the viscosity. 


the vicinity of the shearing instability behaves as 22] 


(e) - (e) 


H 


(e) 


= A e {L-L c )~\ 


(16) 


H 


where (e)n is the value of the average energy far away from the instability, and A e does not 
depend on L. This behavior is clearly observed in the MD simulations, and a fit of the data 
provides the value of L c . In Fig. CD the results obtained for two values of the restitution 
coefficient, a = 0.99 and a = 0.98, and four densities in the range 0.231 < no 2 < 0.385 
are shown. The symbols are from MD simulations and the solid lines are the theoretical 
predictions given by Eq. ©. using the Enskog values for the cooling rate and the shear 
viscosity. It is seen that the Enskog theory provides a good estimation of the critical size. 
A relevant conclusion of this is that the shear viscosity does not have a singular behavior 
when approaching the shearing instability, in agreement with the assumption made in the 
theory developed above. The values of L c obtained in the way described will be used in the 
following to determine the behaviour of the self-difussion coefficient near L = L c . 

The diffusion coefficient can be measured in MD simulations by at least three different 
methods. A first possibility is based on the Einstein relation for the mean square displace- 
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ment of a tagged particle, that for d = 2 is [1 

1 


D = lim —(fr(s) — r 

s—>oo 45 




( 17 ) 


Ail alternative is to use the Green-Kubo formula, Eq. (J5]). The third way of evaluating D 
consists in considering the steady self-diffusion state reached by a mixture of two kinds of 
mechanically equivalent particles, differing only in some label or color, which is in contact 
with two reservoirs for the two types of particles. Then, the diffusion coefficient is obtained 
from the values of the particles flux and the concentration gradient, i.e. from Eq. (j3J) jl7j |. 
The implementations of the three methods have been discussed in detail in the literature 
jl5, 171, and should lead to the same value of the coefficient of self-diffusion if the system 
exhibits diffusive behavior. This consistency was confirmed in our simulations, although 
the results based on the VACF seem to be somewhat more accurate and less noisy, due 
to technical reasons, as it has been already pointed out 16]. For this reason, the values 
being reported in the following were obtained with the Green-Kubo formula. Nevertheless, 
a potential difficulty arises with the tails of the VACF near the shearing instability. For 
large times, the decay is dominated by the hydrodynamic term with the shortest wave 
vector compatible with the boundary conditions, as described by Eq. CD- Given that the 
amplitude of this contribution diverges, it could happen that the exponential long time tail 
would play a significant role in determining the self-diffusion coefficient. In Fig. [2] the long 
time behavior of the VACF is illustrated for a system with a = 0.98, n = 0.3cr -2 , and two 
different system sizes. The exponential decay is clearly identified. Then what has been done, 
instead of using some cut-off in the numerical integration of the VACF to get the diffusion 
coefficient, is to chose for each system a time se > so, such that the exponential decay is 
already clearly observed for that time. For s < se, the results of the simulations have been 
integrated numerically, while for se < s < oo what has been analytically integrated is the 
fit of the simulation data to an exponential, extended up to s —> oo. 

In Fig. [3] the diffusion coefficient measured in a system with n = 0.33ct~ 2 and different 
values of the coefficient of restitution is shown as function of the size of the system. In the 
inelastic cases (a < 1), the divergence shows up when the size of the system approaches 
L c . Far from the divergence, D grows logarithmically with the size. This is the qualitative 
behaviour expected in the elastic limit for all L , as a consequence of the long time tails 
3 J]. Moreover, far from the instability the dependence of the self-diffusion coefficient on the 
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FIG. 2: (Color online) Time evolution of the dimensionless VACF, / = mC(s)/T s t for a system 
of inelastic hard disks with a = 0.98 and n = 0.3a~ 2 . The symbols are simulation results for two 
different system sizes, as indicated in the inset. The time so characterizing the long time regime 
of Eq. (11011 is also given, The solid lines are linear fits identifying the exponential decay for long 
times. 

inelasticity is rather weak. Similar results have been obtained for the other densities studied. 

Also, it has been investigated whether the observed divergence is compatible with the 
logarithmic prediction in Eq. (1151) . Then, in Fig. HJ D is plotted as a function of [— ln(l — 
L/L c )] 1//2 for a system with n = 0.144<r~ 2 and a = 0.95. The linear behavior predicted 
by the theory is observed as the shearing instability is approached, although on a rather 
limited scale range, due to the limitations imposed by the instability of the HCS itself. It 
is interesting to analyze the influence of the mode with the minimum wave vector k m on 
the divergence of D. In the same figure, the contributions to D from the initial decay of 
the VACF up to Sg (A) and from the long time tail described by an exponential (J 2 ), as 
discussed above, are shown. While the former increases quite fast upon approaching L c , 
the latter remains rather small. The reason is that although the amplitude of the time 
exponential associated to k m diverges, as indicated in Eq. CD, its characteristic decay time 
goes to zero, in such a way that the time integral actually remains bounded in the limit 
L —> L c . This has been consistently observed in all the cases investigated. 

In summary, using mode coupling theory the self-diffusion coefficient of a finite granular 
gas in the homogeneous cooling state near the shearing instability has been investigated. 
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FIG. 3: (Color online) Dimensionless diffusion constant D as a function of the size of the system 
L for a = 0.98 (red squares), a = 0.99 (blue circles), and a = 1 (black diamonds). In all cases, it 
is n = 0.33a~ 2 . The critical sizes are lnL c /ir ~ 3.99 for a = 0.98 and lnT c /ir ~ 4.36 for a = 0.99. 
The dashed line is a linear fit of the elastic data. 

In the study, it turns out to be essential to use the macroscopic transport coefficients and 
not the bare (e.g. Enskog) ones. The result given by Eq. (fl4l) implies to substitute D 
by D hyd on the right hand side of Eq. (fT3lh i.e. the transport coefficient considered is the 
(divergent) macroscopic one and not its Enskog value. In this way, it has been predicted that 
the diffusion constant exhibits a logarithmic divergence, that is consistent with the results 
obtained by Molecular Dynamics simulations. The divergence is not dominated by the mode 
with the shortest wavevector, but by a combination of modes in the hydrodynamic region. 
As a by-product, it has been seen that the shear viscosity remains finite at the shearing 
instability. This follows from the fact that this coefficient appears in the hydrodynamic 
expression of the critical size for the system to exhibit the unstable behaviour. If the viscosity 
would be infinite, the instability never would be observed in practice. The same property 
is crucial for the derivation of Eq. (TT5l) from Eq. (TT3lh and the latter has been shown to be 
consistent with the behaviour observed in MD simulations. More extensive results will be 
published elsewhere. It is worth to emphasize that the mode coupling calculation used here, 
relies solely on the hydrodynamic modes and, therefore, its confirmation by the simulation 
results is another piece of evidence for the relevance and importance of hydrodynamics in 
these systems, a fact often discussed in the literature. 
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FIG. 4: (Color online) Dimensionless diffusion constant D = I± + Io (black stars) as a function 
of (—InA) 1 / 2 , where A = (L c — L)/L c , for a system with a = 0.95 and n = 0.144ct - 2 . The solid 
line is a linear fit corresponding to the theoretical prediction, Eq. (|15l) . The component Jo (red 
squares) is the contribution from the long time exponential decay of the VACF, while I\ (blue 
circles) is the contribution from shorter times. 

This research was supported by the Ministerio de Economia y Competitividad (Spain) 
through Grant No. FIS2014-53808-P (partially financed by FEDER funds). 


[1] H.M. Jaeger, S.R. Nagel, and R.P. Behringer, Rev. Mod. Phys. 68, 1259 (1996). 

[2] I. Goldhirsch, Annu. Rev. Fluid Mech. 35, 267 (2003). 

[3] B.J. Alder and T.E. Wainwright, Phys. Rev. Lett. 18, 988 (1967). 

[4] J.R. Dorfman, in Fundamental Problems in Statistical Mechanics, Vol. 3, E.D.G. Cohen, editor 
(North-Holland Publishing Company, Amsterdam, 1975). 

[5] P. K. Haff, J. Fluid Mech. 134, 401 (1983). 

[6] J.J. Brey, J.W. Dufty, C.S. Kim, and A. Santos, Phys. Rev E 58, 4638 (1998). 

[7] V. Garzo and J.W. Dufty, Phys. Rev. E 59, 5895 (1999). 

[8] A. Baskaran, J.W. Dufty, and J.J. Brey, Phys. Rev. E 77, 031311 (2008). 

[9] I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993). 

[10] S. McNamara and W.R. Young, Phys. Rev. E 50, R28 (1994). 


10 








[ 11 ] 

[ 12 ] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[ 20 ] 
[ 21 ] 
[ 22 ] 

[23] 

[24] 

[25] 


J.J. Brey, M.J. Ruiz-Montero, and D. Cubero, Phys. Rev. E 60, 3150 (1999). 

J.J. Brey, M.J. Ruiz-Montero, and A. Dominguez, Phys. Rev. E 78, 041301 (2008). 

J.J. Brey, M.J. Ruiz-Montero, D. Cubero, and R. Garcia Rojo, Phys. Fluids 12, 876 (2000). 
N. V. Brilliantov and T. Poschel, Phys. Rev. E 61, 1716 (2000). 

J.W. Dufty, J.J. Brey, and J. Lutsko, Phys. Rev. E 65, 051303 (2002). 

J. Lutsko, J.J. Brey, and J.W. Dufty, Phys. Rev. E 65, 051304 (2002). 

J.J. Brey and M. J. Ruiz-Montero, Phys. of Fluids 25, 113302 (2013). 

J.F. Lutsko, Phys. Rev. E 63, 061211 (2001). 

J. J. Brey, M.J. Ruiz-Montero, and F. Moreno, Phys. Rev. E 69, 051303 (2004). 

J. J. Brey and M.J. Ruiz-Montero, Phys. Rev. E 91, 012202 (2015). 

V. Garzo, A. Santos, and J.M. Montanero, Physica A 376, 94 (2007). 

J. J. Brey, A. Dominguez, M. I. Garcia de Soria, and P. Maynar, Phys. Rev. Lett. 96, 158002 


(2006). 

A factor of 2 is missing in Eq. (68) of Ref. 


M- 


H. Hayakawa and M. Otsuki, Phys. Rev. E 76, 051304 (2007). 

A. Fiege, T. Aspelmeier, and A. Zippelius, Phys. Rev. Lett. 102, 098001 (2009). 


11 



